	program cst
      include 'cst.inc'      
c****************************************************************************************
c	numerical model of large-scale coastal dynamics based on the
c	coastal system tract concept (see cowell et al 2000). this is
c	a community-developed hybrid model on time-average coastal and
c     processes.
c****************************************************************************************
      character cht*4
c
c --- Read parameters ---
c
      call read_para
c
c --- Initialize the cstat and d terms ---
c
	call cd_init
c
c --- Save initial data for postprocessing ---
c     nn is for time step
c
	nn=1
	do i=1,imax
	   syser(i,nn)=yshore(i)
	   sxser(i,nn)=((float(i))*deltax)-(deltax/2.)
     	   do j=1,jmax
	      hseries(i,j,nn)=h(i,j)
	   enddo
	enddo      
c
c --- set the parameters for subroutine shoreline ---
c     hold is updated in the shoreline subroutine
c
	do i= 1,imax
         hold(i)=h(i,jshore(i))
      enddo   
c
	deltah=htop-hbot
	tanalpha=deltah/deltay
	ta=tanalpha/2.
	hav=htop-(deltah/2.)
	ta=tanalpha/2.
	haa=((-hbot)**2)/(2.*deltah)
	hbb=(htop**2)/(2.*deltah)
	ha=hbot+haa
	hb=htop-hbb
c
	ylamax=abs(hbot/tanalpha)
	ylbmax=abs(htop/tanalpha)
c
	fma=abs((ylamax)/deltay)
	fmb=abs((ylbmax)/deltay)
c***************************************
c*    begin the main time loop         *
c***************************************
	iisave=0
	totalin=0.0  !Added by Ming Ye on 7/25/2009
	totalout=0.0 !Added by Ming Ye on 7/25/2009
	sumqin=0.0   !Added by Ming Ye on 7/25/2009
	sumqout=0.0  !Added by Ming Ye on 7/25/2009
	sumbot=0.0   !Added by Ming Ye on 7/25/2009
	do m=1,itmax
   	   write(*,*)'itmax= ',m
	   iisave=iisave+1
c
c --- Update parameters of the large-scale diffusion, Equation (4) of Niedoroda (1995) ---
c
	   call c_d
c	   write(cht,'(i4.4)')m
c	   open(35,file='c_d_'//cht//'.dat')
c	   do i=1,imax
c	      do j=1,jmax
c	        write(35,'(3i5,3f9.2)')i,j,jshore(i),dl(i,j),dc(i,j),cstar(i,j)
c	      enddo
c	   enddo
c	   close(35)        
c
c --- Calculate the angle between the tangent line and negative x direction --- 
c
	   call newangle
c
c --- Update Wave propagation and refraction near shore ---
c
	   call wave_prop
c
c --- Calculate sediment transport in the surf zone ---- 
c
	   call surfzone(m)
c
c --- Calculate sedoment transport in the shelf, the Qi term, Equation (2) of Niedoroda (1995) ---
c
	   call shelf_qs
c
c --- Solve for the exner continuity equation of depth h, Equations (10) and (11) of Niedoroda (1995) --- 
c
         call exner
c
c --- Solve for the new shoreline position ---
c
	   call shoreln4  
c
	   do i=1,imax
	      yshore(i)=yoff(i)+(jshore(i)-0.5)*deltay
         enddo
c
c --- Track sediment flow in and out of system ---
c
 	   totalin=totalin+qlit(1)*deltay*dt*(1./(1.-pore))
	   totalout=totalout+qlit(imax-1)*deltay*dt*(1./(1.-pore))
c
c --- Calculate 
c 	   
	   sumqint=0.0
	   sumqoutt=0.0
	   do j=1,jmax
	      sumqint=sumqint+qxrt(1,j)*deltay*dt*(1./(1.-pore))
	      sumqoutt=sumqoutt+qxrt(imax-1,j)*deltay*dt*(1./(1.-pore))
         enddo
	   sumqin=sumqin+sumqint
	   sumqout=sumqout+sumqoutt
	   sumbott=0.0
	   do i=1,imax
	      sumbott=sumbott+qyup(i,1)*deltax*dt*(1./(1.-pore))
         enddo
	   sumbot=sumbot+sumbott

	   if(iisave.eq.nnsave) then 
	     iisave=0
	     nn=nn+1
	     ndata=nn
	     do i=1,imax
	        syser(i,nn)=yshore(i)
	        sxser(i,nn)=((float(i))*deltax)-(deltax/2.)
     	        do j=1,jmax
	           hseries(i,j,nn)=h(i,j)
	        enddo
	     enddo
	   endif        
      enddo
c
c	calulate mass change
      delh=0  !Added by Ming Ye on 7/25/2009
      delhmax=-1.0E+20  !Added by Ming Ye on 7/25/2009
      delhmin=1.0E+20   !Added by Ming Ye on 7/25/2009
	do i=2,imax-1
	   do j=2,jmax
	      delht=h(i,j)-hseries(i,j,1) 
	      delhmax=amax1(delhmax,delht)
	      delhmin=amin1(delhmin,delht)
	      delh=delh+delht
	   enddo   
      enddo
      
	write(*,*)'mass balance:'
	write(*,*)'lit from left: = ',totalin
	write(*,*)'lit off right: = ',totalout
	write(*,*)'qin from left: = ',sumqin
	write(*,*)'qout off right: = ',sumqout
	write(*,*)'qin from bot : = ',sumbot
	write(*,*)'accum from lit = ',totalin-totalout
	write(*,*)'accumlation = ',delh*deltax*deltay
	write(*,*)'max change in height = ',delhmax
	write(*,*)'min change in height = ',delhmin

c
	open(unit=2,file='contrl.dat')
	write(2,*)imax,jmax,ndata,deltax,deltay,dt
	close(unit=2)
c
	open(unit=2,file='cstmout.dat')
      do n=1,ndata
	   do j=1,jmax
   	      write(2,*)(hseries(i,j,n),i=1,imax)
   	   enddo   
      enddo
	do n=1,ndata
   	   write(2,*)(sxser(i,n),syser(i,n),i=1,imax)
   	enddo   
	close(unit=2)
c
c
c	the following is for diagnoses only
c
	open(unit=2,file='dgrid_g.dat')
	write(2,*)'variables = "x", "y", "z", "theta"'
	write(2,*)'zone f=point, i=',imax,' ,j=',jmax
	do j=1,jmax
	   do i=1,imax
	      dddx=deltax/2.
	      dddy=deltay/2.
 	      write(2,*)(deltax*(i-1))+dddx,(deltay*(j-1))+dddy,h(i,j),a(i,j)
 	   enddo
 	enddo      
	close(unit=2)

 	open(unit=2,file='hgtchng.dat')
	write(2,*)'variables = "x", "y", "z", "theta"'
	write(2,*)'zone f=point, i=',imax,' ,j=',jmax
	do j=1,jmax
	   do i=1,imax
	      dddx=deltax/2.
	      dddy=deltay/2.
 	      write(2,*)(deltax*(i-1))+dddx,(deltay*(j-1))+dddy,
     +      h(i,j)-hseries(i,j,1),a(i,j)    
         enddo
      enddo    
      close(2)

      stop	
      end